#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <algorithm>

#define RUNS 10

static void HandleError(cudaError_t err,
                        const char *file,
                        int line) {
    if (err != cudaSuccess) {
        printf("%s in %s at line %d\n", cudaGetErrorString(err),
               file, line);
        exit(EXIT_FAILURE);
    }
}

#define HANDLE_ERROR(err) (HandleError(err, __FILE__, __LINE__))
#define HANDLE_NULL(a) {if (a == NULL) { \
                        printf("Host memory failed in %s at line %d\n", \
                               __FILE__, __LINE__); \
                        exit(EXIT_FAILURE);}}

const int N = 1024 * 1024 * 64;     //积分时划分的份数
const int threadsPerBlock = 256;    //block中的线程数
const int blocksPerGrid = 64;       //grid中的block数

//缺省__host__，表明CPU运行，CPU调用
double function_for_cpu(double x) {
    return 4 / (1 + x * x);
}

//__device__修饰，只能被内核函数调用，表明GPU运行，GPU调用
__device__ double function_for_gpu(double x) {
    return 4 / (1 + x * x);
}

/**
 * __global__修饰，内核函数，表明GPU运行，CPU调用
 * 在GPU上使用积分法并行计算PI
 * @param a 积分下界
 * @param b 积分上界
 * @param integral 存储积分结果
 * @return
 */
__global__ void trap(double *a, double *b, double *integral) {

    /**
     * __shared__修饰，声明一个共享内存缓冲区，名字为cache
     * 表示数据存放在共享存储器中，每个线程块都有该变量的一个副本；
     * 只有在块内的线程可以访问，其它块内的线程不能访问；
     */
    __shared__ double cache[threadsPerBlock];

    /**
     * 该步的目的是计算初始线程索引，其中threadIdx, blockIdx, blockDim都是内置变量
     * threadIdx是存储线程信息的结构体，对于线程0来说，threadIdx.x=0;
     * blockDim.x表示block在x维度的线程数量，本例中使用的是一维线程块，
     *     因此只需用到blockDim.x
     * blockIdx.x表示block的索引，对于第一个线程块来说，blockIdx.x=0;
     *     对于第二个线程块来说，blockIdx.x=1...
     * 在计算tid线程索引时，需要要在threadIdx.x 的基础上加上一个基地址，
     *     实际上就是将二维索引空间转换为线性空间
     */
    int tid = threadIdx.x + blockIdx.x * blockDim.x;

    /**
     * 共享内存缓存中的偏移就等于线程索引，线程块索引与该偏移无关，
     * 因为每个线程块都拥有该共享内存的私有副本
     */
    int cacheIndex = threadIdx.x;

    double x, temp = 0;
    while (tid < N) {
        x = *a + (double)(*b - *a) / N * (tid + 0.5);
        temp += function_for_gpu(x);

        /**
         * 在每个线程计算完当前索引上的任务后，需要对索引进行递增，
         * 其中，递增的步长为线程格中正在运行的线程数量，
         * 这个数值等于线程块中的线程数量乘以线程格中线程块的数量，
         * 即 blockDim.x * gridDim.x
         *
         * 该方法类似于多CPU或多核CPU的并行，数据迭代的增量不是1，
         *     而是CPU的数量；在GPU实现中，一般将并行线程数量看做处理器的数量
         */
        tid += blockDim.x * gridDim.x;
    }

    //设置cache中相应位置上的值
    cache[cacheIndex] = temp;

    /**
     * 对线程块中的线程进行同步，该操作用于确保
     *     所有对共享数组cache[]写入操作在读取cache之前完成
     */
    __syncthreads();

    /**
     * 对于归约运算来说，以下代码要求threadsPerBlock必须是2的幂，
     * 因为每次合并，要求分成的两部分数组的长度要一致
     * 基本思想：
     *     每个线程将cache[]中的两个值相加起来，然后将结果保存回cache[]
     *     由于每个线程都将两个值合并为一个值，那么在完成这个步骤后，
     *     得到的结果数量就是计算开始时数值数量的一半。接着，对这一半
     *     进行相同操作，直到cache[]中256个值归约为1个值
     */
    int i = blockDim.x / 2;
    while (i != 0) {
        if (cacheIndex < i)
            cache[cacheIndex] += cache[cacheIndex + i];
        __syncthreads();
        i /= 2;
    }

    //将这个值保存到全局内存后，内核函数结束
    //这里使用了索引为0的线程将cache[0]写入全局内存
    if (cacheIndex == 0)
        integral[blockIdx.x] = cache[0];
}

// 显存占用查询函数
void printGpuMemoryUsage(const char* label) {
    size_t free_mem, total_mem;
    cudaMemGetInfo(&free_mem, &total_mem);
    printf("[%-15s] GPU Memory Usage: Used = %.2f MB, Free = %.2f MB, Total = %.2f MB\n",
           label,
           (total_mem - free_mem) / (1024.0 * 1024.0),
           free_mem / (1024.0 * 1024.0),
           total_mem / (1024.0 * 1024.0));
}

int main(void) {
    double integral;
    double a = 0, b = 1;
    size_t total_mem, free_mem_before, free_mem_after;
    size_t peak_mem_global = 0;
    float total_time = 0.0f;

    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    for (int i = 0; i < RUNS; ++i) {
        double *partial_integral = (double*)malloc(blocksPerGrid * sizeof(double));
        double *dev_partial_integral;
        double *dev_a, *dev_b;

        // 显存峰值变量，单次峰值
        size_t peak_mem = 0;

        // 记录显存分配前
        cudaMemGetInfo(&free_mem_before, &total_mem);

        // 申请GPU内存
        HANDLE_ERROR(cudaMalloc((void**)&dev_a, sizeof(double)));
        HANDLE_ERROR(cudaMalloc((void**)&dev_b, sizeof(double)));
        HANDLE_ERROR(cudaMalloc((void**)&dev_partial_integral, blocksPerGrid * sizeof(double)));

        // 显存分配后更新峰值
        cudaMemGetInfo(&free_mem_after, &total_mem);
        peak_mem = std::max(peak_mem, total_mem - free_mem_after);

        // 复制数据到GPU
        HANDLE_ERROR(cudaMemcpy(dev_a, &a, sizeof(double), cudaMemcpyHostToDevice));
        HANDLE_ERROR(cudaMemcpy(dev_b, &b, sizeof(double), cudaMemcpyHostToDevice));

        // 内核调用前显存检测
        size_t mem_before_kernel;
        cudaMemGetInfo(&mem_before_kernel, &total_mem);
        peak_mem = std::max(peak_mem, total_mem - mem_before_kernel);

        // 计时开始
        cudaEventRecord(start, 0);

        trap<<<blocksPerGrid, threadsPerBlock>>>(dev_a, dev_b, dev_partial_integral);

        // 计时结束
        cudaEventRecord(stop, 0);
        cudaEventSynchronize(stop);

        // 内核调用后显存检测
        size_t mem_after_kernel;
        cudaMemGetInfo(&mem_after_kernel, &total_mem);
        peak_mem = std::max(peak_mem, total_mem - mem_after_kernel);

        // 拷贝回CPU
        HANDLE_ERROR(cudaMemcpy(partial_integral, dev_partial_integral, blocksPerGrid * sizeof(double), cudaMemcpyDeviceToHost));

        // CPU归约
        integral = 0;
        for (int j = 0; j < blocksPerGrid; ++j) {
            integral += partial_integral[j];
        }
        integral *= (b - a) / N;

        float ms;
        cudaEventElapsedTime(&ms, start, stop);

        // 释放显存
        HANDLE_ERROR(cudaFree(dev_a));
        HANDLE_ERROR(cudaFree(dev_b));
        HANDLE_ERROR(cudaFree(dev_partial_integral));
        free(partial_integral);

        // 记录峰值显存（释放后检查）
        size_t free_mem_end;
        cudaMemGetInfo(&free_mem_end, &total_mem);
        peak_mem = std::max(peak_mem, total_mem - free_mem_end);
        peak_mem_global = std::max(peak_mem_global, peak_mem);

        printf("Run %3d: GPU time = %.6f ms, Peak memory usage = %.2f MB, Integral = %.10lf\n", i, ms, peak_mem / (1024.0 * 1024.0), integral);

        // 累计第1~99次时间
        if (i > 0) {
            total_time += ms;
        }
    }

    printf("\nAverage GPU time over %d runs (excluding 1st run): %.6f ms\n", RUNS - 1, total_time / (RUNS - 1));
    printf("Global peak GPU memory usage over all runs: %.2f MB\n", peak_mem_global / (1024.0 * 1024.0));

    cudaEventDestroy(start);
    cudaEventDestroy(stop);

    return 0;
}
